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ABSTRACT 

In this paper, we consider the dynamics of a system of hot super-Earths 
or Neptunes such as HD 40307. We show that, as tidal interaction with the 
central star leads to small eccentricities, the planets in this system could be 
undergoing resonant coupling even though the period ratios depart signifi- 
cantly from very precise commcnsurability. In a three planet system, this is 
indicated by the fact that resonant angles librate or are associated with long 
term changes to the orbital elements. In HD 40307, we expect that three res- 
onant angles could be involved in this way. We propose that the planets in 
this system were in a strict Laplace resonance while they migrated through 
the disc. After entering the disc inner cavity, tidal interaction would cause the 
period ratios to increase from two but with the inner pair deviating less than 
the outer pair, counter to what occurs in HD 40307. 

However, the relationship between these pairs that occurs in HD 40307 
might be produced if the resonance is impulsively modified by an event like 
a close encounter shortly after the planetary system decouples from the disc. 
Wc find this to be in principle possible for a small relative perturbation on the 
order of a few x 10~ 3 but then we find that evolution to the present system in 
a reasonable time is possible only if the masses are significantly larger than the 
minimum masses and tidal dissipation is very effective. On the other hand we 
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found that a system like HD 40307 with minimum masses and more realistic 
tidal dissipation could be produced if the eccentricity of the outermost planet 
was impulsively increased to ~ 0.15. 

We remark that the form of resonantly coupled tidal evolution we consider 
here is quite general and could be of greater significance for systems with inner 
planets on significantly shorter orbital periods characteristic of for example 
CoRoT 7 b. 

Key words: planetary systems: formation — planetary systems: protoplan- 
etary discs 

1 INTRODUCTION 

Neptune mass extrasolar planets around main sequence stars were first detected five years 
ago. Since then, 27 planets with a projected mass lower than 25 earth masses have been 
discovered, the lightest one having a projected mass of 1.94 M®. Among these objects, 17 
have a semi-major axis smaller than 0.1 astronomical unit (au). They are called hot Neptunes 
or hot super-Earths. So far, two multiple planet systems with at least two such objects have 
been observed. One of them is a four planet system around the M-dwarf GJ 581 (Bonfils 
et al. 2005, Udry et al. 2007, Mayor et al. 2009a). The projected masses of the planets are 
1.9, 15.6, 5.4 and 7.1 M® and the periods are 3.15, 5.37, 12.93 and 66.8 days, respectively. 
Note that although the eccentricities of the two innermost planets are compatible with 
zero, those of the next outermost and outermost planets are 0.16 and 0.18 respectively. The 
second multiple system which we shall focus upon in this paper is that around the K-dwarf 
HD 40307 (Mayor et al. 2009b). It comprises three planets with projected masses of 4.2, 6.9 
and 9.1 M® and periods of 4.31, 9.62 and 20.46 days, respectively. The eccentricities are all 
compatible with zero. 

If we label the planets in a system with successive integers starting from the innermost 
labelled, 1, and moving outwards, then, in the system Gl 581, the ratio of the periods of 
planets 2 and 1 is 2.41 and that of the periods of planets 3 and 2 is 1.70. We note that 
these numbers depart from 5/2 and 5/3 only by 4% and 2%, respectively. a In the HD 40307 
system , the ratio of the periods of planets 2 and 1 is 2.23 and that of planets 3 and 2 is 2.13, 
which depart from 2 by 11.5% and 6.5%, respectively. Because departure from exact mean 
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motion resonances in these systems is significant, the importance of resonances for their 
dynamical evolution has been ruled out (Mayor et al. 2009a, 2009b, Barnes et al. 2009). We 
shall address this aspect for the HD 40307 system in this paper. 

Migration due to tidal interaction with the disc is probably the mechanism by which 
planets end up on short period orbits, as in situ formation requires very massive discs (e.g., 
Raymond et al. 2008). Planet-planet scatterings may also lead to short period orbits, but 
only when tidal circularization is efficient. According to Kennedy & Kenyon (2008), scatter- 
ings are not a likely way of producing hot super-Earths, because of the long circularization 
timescales involved. In a previous paper (Terquem & Papaloizou 2007, see also Brunini & 
Cionco 2005), we proposed a scenario for forming hot super-Earths on near commensurable 
orbits, in which a population of cores that assembled at some distance from the central star 
migrated inwards due to tidal interaction with the disc while merging on their way in. We 
found that evolution of an ensemble of cores in a disc almost always led to a system of a few 
planets, with masses that depended on the total initial mass of the system, on short period 
orbits with mean motions that frequently exhibited near commensurabilities and, for long 
enough tidal circularization times, apsidal lines that were locked together. Starting with a 
population of 10 to 25 planets of 0.1 or 1 M e , we ended up with typically between two and 
five planets with masses of a few tenths of an earth mass or a few earth masses inside the 
disc inner edge. Interaction with the central star led to the initiation of tidal circulariza- 
tion of the orbits which, together with possible close scatterings and final mergers, tended 
to disrupt mean-motion resonances that were established during the migration phase. The 
system, however, often remained in a configuration in which the orbital periods were close to 
commensurability. Apsidal locking of the orbits, if established during migration, was often 
maintained through the action of these processes. 

This scenario has been questioned (Mayor et al. 2009a, 2009b) because the multiple 
planet systems of hot super-Earths detected so far do not exhibit close enough mean mo- 
tion commensurabilities. In this paper, we argue that the system around HD 40307 does 
actually exhibit the effects of resonances through secular effects produced by the action of 
the resonant angles coupled with the action of tides from the central star, and that such a 
configuration could be produced if the system formed as described in Terquem & Papaloizou 
(2007), but with the addition of relatively small perturbations to the orbits arising from close 
encounters or collisions after the system enters the disc inner cavity. The reason that what 
are regarded as resonant effects can arise even when departures from strict commensurability 
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are apparently large is that tidal circularization produces small eccentricities which, for first 
order resonances, can be consistent with resonant angle libration under those conditions (see 
Murray & Dermott 1999). 

We consider in detail the evolution of this system and other putative similar systems with 
scaled masses and orbital periods after they form a strict three planet Laplace resonance 
as a result of convergent inward orbital migration. We go on to consider the onset of tidal 
circularization and how it leads to a separation of the semi- major axes and consequent 
increasing departure from commensurability. We find that if the strict Laplace resonance is 
broken by a fairly small relative perturbation corresponding to a few parts in a thousand, 
continuing evolution in a modified form of the three planet resonance could in principle in the 
case of HD 40307 lead to a system like the observed one. However, this would require masses 
significantly larger than the minimum estimates and possibly unrealistically efficient tidal 
dissipation with the tidal dissipation parameter Q' ~ 10. However, weaker tidal effects could 
be responsible for some of the deviation from strict commensurability, with the rest being 
produced by a larger perturbation resulting from the processes mentioned above inducing 
eccentricities of the order 0.1. Note however that tidal evolution could play a more significant 
role in similar systems with shorter orbital periods. 

In section |2j we describe the numerical model we use to simulate the evolution of a 
system of three planets migrating in a disc and evolving under the tidal interaction with 
the star after they enter the disc inner cavity. In section [3j we simulate the system around 
HD 40307 using either the minimum or twice the minimum masses for the planets under 
varying assumptions about initial eccentricities and the strength of the tidal interaction. We 
show that three of the four possible resonant angles associated with 2:1 commensurabilities 
either librate or have long term time averages, indicating evolution of the system towards a 
resonant state driven by the tidal influence of the star. 

In a related appendix, we give a semi-analytic model for a system in a three planet 
resonance undergoing circularization with departures from strict mean motion commensu- 
rabilities. We note in passing that the discussion may also be applied to a two planet system 
near a 2:1 resonance. We show that the departure from commensurability increases with time 
being oc t 1 / 3 and derive an expression for the timescale required to attain a given departure 
that can be used to interpret extrapolate from and generalize the numerical simulations. 

In section |4j we present results of numerical simulations of a three planet system migrat- 
ing in a disc, establishing a 4:2:1 Laplace like resonance, and evolving under tidal effects 
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induced by the star after entering the disc inner cavity. As expected, as the orbits are being 
tidally circularized, the planets move away from exact commensurabilites, and the period 
ratios of the outer and inner pairs of planets increases. 

The system is still resonant in the sense that some of the resonant angles continue to 
librate. To get a larger period ratio for the inner pair of planets, as in HD 40307, some 
disruption of this Laplace resonant state is needed. In our simulations, this is achieved by 
applying an impulse (that could result from an encounter or collision) to the system. This 
leads to departures from exact commensurabilites of the form seen in HD 40307 while still 
retaining libration of three of the four resonant angles. Finally, in section [5] we discuss and 
summarize our results. 



2 MODEL AND INITIAL CONDITIONS 

The model we use has been described in Terquem & Papaloizou (2007). We recall here the 
main features. We consider a system consisting of a primary star and N planets embedded 
in a gaseous disc surrounding it (for HD 40307, N =3). The planets undergo gravitational 
interaction with each other and the star and are acted on by tidal forces from the disc 
and the star. The evolution of the system is studied by means of the solution of iV-body 
problems, in which the tidal interactions are included as dissipative forces. 
The equations of motion are: 

(Prs GM*Yi Gin,-, (v: — r,-) „ „ „ 

■7«F = --rir- E i |3 -r + r t + r r , 1 

at 1 r,- d ' r,- — r.-r 

where M*, rrij and r 3 - denote the mass of the central star, that of planet j and the position 
vector of planet j, respectively. The acceleration of the coordinate system based on the 
central star (indirect term) is: 

j=l I 3\ 

and that due to tidal interaction with the disc and/or the star is dealt with through the 

addition of extra forces as in Papaloizou & Larwood (2000): 

1 dti 2 / dYi \ 2 / d,Ti 



t mji dt \r t \H e , \dt 'V" 1 hi \dt '*')*" (3) 
where t m ^, t e ^ and t^i are the timescales over which, respectively, the angular momentum, 

the eccentricity and the inclination with respect to the unit normal e 2 to the assumed fixed 

gas disc midplane change. Evolution of the angular momentum and inclination is due to 
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tidal interaction with the disc, whereas evolution of the eccentricity occurs due to both tidal 
interaction with the disc and the star. We have: 

T- = TT + tt, (4) 

L e,% L e ,i L e,i 

where t d e i and t s e i are the contribution from the disc and tides raised by the star, respectively. 
Relativistic effects are included through T r ( see Papaloizou & Terquem 2001). 

The way it is implemented, interaction with the star does not modify the angular mo- 
mentum of a single orbit (the orbital decay timescale due to tidal interaction with the star, 
which as indicated below is estimated to be much longer than the circularization timescale, 
has been ignored here). However, in the formulation above, eccentricity damping causes ra- 
dial velocity damping, which results in energy loss at constant angular momentum. As a 
consequence, both the semi-major axis and eccentricity are reduced together (at least for a 
single planet). 

2.1 Orbital circularization due to tides from the central star 

The circularization timescale due to tidal interaction with the star can be estimated, from 
expressions given by Goldreich & Soter (1966), to be: 

M ffi \ 2/3 /20a, A 6 ' 5 



t s e>i = 4.065 x 10 4 — M Q' years, (5) 

V rrii J \1 an J 

where cij is the semi-major axis of planet i. Here we have adopted a mass density of 1 g cm 3 
for the planets (uncertainties in this quantity could be incorporated into a redefinition of 
Q'). The parameter Q' = 3Q/(2/c2), where Q is the tidal dissipation function and k 2 is the 
Love number. Equation (|5| applies to the situation where tides raised by the central star are 
dissipated in the planet in the limit of small eccentricity which is appropriate here. Tides 
raised by the planet in the star have been estimated to be unimportant for the low mass 
planets of interest (eg. Barnes et al 2009). For solar system planets in the terrestrial mass 
range, Goldreich & Soter (1966) give estimates for Q in the range 10-500 and k 2 ~ 0.3, 
which correspond to Q' in the range 50-2500. Of course this parameter must be regarded as 
being very uncertain for extrasolar planets, and accordingly we have considered values of Q' 
in the range 10 — 1000 in this paper. For an earth mass planet at 0.05 au and Q' = 50, we 
get t s ei = 2 x 10 6 years. This indicates that circularization processes will be important for 
close in extrasolar planets in this mass range over the entire range of Q' considered. 

But note that although circularization operates on each orbit according to equation (|5]), 
because the planets interact they may enter into a decay mode with a decay timescale that is 
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not straightforwardly estimated as above (see section 3.1) . Nonetheless the above estimate 
is found to be indicative. 

We have pointed out above that, in our formulation, tidal interaction with the star does 
not change the angular momentum of a single orbit but modifies the semi-major axis. The 
physical basis for this is that the planets rapidly attain pseudo-synchronization (e.g., Ivanov 
& Papaloizou 2007), after which they cannot store significant angular momentum changes 
through modifying their intrinsic angular momenta. 

We restrict consideration for the moment to a single planet. Since the orbital angular 
momentum is conserved, (ddi / 'dt) / ' a\ ~ —2ef/t s ei for small eccentricities, where e; is the 
eccentricity of the planet we consider. Therefore a, changes only by about 5 percent in 
10 10 years for Q' = 10, m; = 1 M® and e\ = 10~ 3 . However, larger changes may occur for 
larger planet masses. We investigate the dependence of this form of evolution on the system 
configuration in detail below. 

Finally we remark that all the simulation results reported here were obtained assuming 
that Q' is the same for all the planets. However, tests were carried out with this assumption 
relaxed and these are alluded to in section |4] below. 

2.2 Type I migration 

In the local linear treatment of type I migration (e.g., Tanaka et al. 2002), if the planet is 
not in contact with the disc, there is no interaction between them so that t m> i, t d ei and 
are taken to be infinite. When the planet is in contact with the disc, disc-planet interactions 
occur leading to orbital migration as well as eccentricity and inclination damping (e.g., 
Ward 1997). In that case, away from the disc edge, for some simulations we have adopted: 

e,- 



146.0 



1 + 



1.3H/r 



K,i = 0.362 



1 + 0.25 ' 



l.lH/r 

3 



H/ry M M e a, 

years, (6) 



H/r 



0.05 J M d rrii 1 au 
H/ r y M M ffi a % 



0.05 / M d rrii I au 



years, (7) 



and tj j = t e> i (eq. [31] and [32] of Papaloizou & Larwood 2000 with f s = 0.6). Here e, is 
the eccentricity of planet i, H/r is the disc aspect ratio and M d is the disc mass contained 
within 5 au. We have assumed here that the disc surface mass density varies like r~ 3 / 2 . For a 
1 earth mass planet on a quasi-circular orbit at 1 au, we get t m ^ ~ 10 5 yr and t d ei ~ 500 yr 
for M d = 10- 3 M and H/r = 0.05. The timescales given by equations ^ and ^ can be 
used not only for small values of e», but also for eccentricities larger than H/r. 
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We have also carried out simulations with t mj i and t ei taken to be constants and thus 
independent of eccentricity. The commensurabilities that are formed in the system depend on 
the ratio of the migration rate to the orbital frequency, with commensurabilities of low order 
and low degree forming when this ratio is small. Therefore, using equations (|6]) and ^ or 
constant timescales give very similar results provided the eccentricity damping is ultimately 
effective and the magnitudes of the migration rates are comparable in both cases at zero 
eccentricity. 

The type I migration rates to be used are uncertain even when the surface density is 
known, largely because of uncertainties regarding the effectiveness of coorbital torques (e.g., 
Paardekooper & Melema 2006, Pardekooper &; Papaloizou 2008, 2009). In this context there 
are indications from modelling the observational data that the adopted type I migration 
rate should be significantly below that predicted by ^ and <Q or Tanaka et al (2002) (see 
Schlaufman et al 2009). Hence we explore a range of possible values. 

3 HD 40307 AND THREE PLANET RESONANCE 

Three super-Earths have been detected orbiting the star HD 40307 through the radial 
velocity variations they produce (Mayor et al. 2009a). The planets have projected masses 
of 4.2, 6.9 and 9.2 and the reported semi-major axes are 0.047, 0.081 and 0.134 au, 
respectively, which correspond to period ratios of 2.26 and 2.13 for the inner and the outer 
pairs of planets, respectively. The periods reported by Mayor et al. (2009a) are 4.31, 9.62 and 
20.46 days, which correspond to period ratios of 2.23 and 2.13. The derived eccentricities are 
compatible with zero. We will refer to the planets as 1, 2 and 3, planet 1 being the innermost 
object. 

Because the period ratios differ significantly from 2, resonances have been ruled out 
(Mayor et al. 2009a, Barnes et al. 2009). However, the period ratios are not far from 4:2:1, 
being characteristic of a three planet resonance. This indicates that such a resonance might 
play a role at some evolutionary phase of the system. To investigate this hypothesis, we have 
simulated the system by means of iV-body calculations using the method described in the 
previous section. The equations of motion are integrated using the Bulirsch-Stoer method 
(e.g., Press et al. 1993). 

The planets were set up on circular orbits around a 0.77 M Q star (the estimated mass of 
HD 40307) with both the minimum masses indicated above and also masses a factor of two 
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larger. The initial semi-major axes were taken to be those indicated above, except in two 
simulations described below, where the innermost semi-major axis was adjusted to make 
the period ratio of the innermost pair 2.23 instead of 2.26. 

The results of a calculation with the tidal parameter Q' — 10 and minimum masses 
are displayed in Fig. [TJ which shows the evolution of the eccentricities and runnimg time 
averages of the cosines of the resonant angles <£>i = 2A 2 — Ai — w\, $ 2 = 2A 2 — Ai — w 2 , 
$3 = 2A3 — A 2 — W2 and $4 = 2A3 — A 2 — -073 over a timescale of 3.6 million years. Here 
Aj and cDj are the mean longitude and the argument of pericentre of planet i, respectively. 
The semi-major axes do not vary significantly with time while the eccentricities remain 
small, below 1CT 3 . However, tidal circularization processes are important for this run as the 
eccentricity of the outermost planet and the amplitude of the oscillations of the eccentricity 
of the inner two planets tends to decrease. Interestingly, the running time averages of the 
cosines of the resonant angles are not zero, as would be expected for a non resonant system. 
They are such that the values for cos $1 and cos $3 approach unity, while the value for cos $4 
moves towards —1, corresponding to the angle approaching 7r. The angle $ 2 exhibits non 
resonant behaviour. 

We have also investigated the dependence of the evolution on Q', small variations of the 
initial semi-major axes and the magnitudes of the planet masses. Figure [2] shows the evolution 
for Q' = 1000, minimum masses and the initial innermost semimajor axis adjusted to make 
the innermost period ratio 2.23. In this case, the mean values of the eccentricities show 
evidence only of limited evolution. Nonetheless, running time averages of the cosines of $1 
$3 and $4 quickly attain steady non zero values. This is consistent with the evolution of these 
runs being driven by tidal effects and towards a resonant state in which the resonant angles 
are associated with secular evolution of the orbital elements. This is able to occur without 
precise commensurability only because the eccentricities are very small (see appendix A). 
For comparison in Fig. [3] we plot the running means of eicos$i, e 2 cos$3, and e3Cos$4 
for the same run. These also attain non zero values also indicating long term variations 
associated with the corresponding resonant angles. Note that the cosines of the resonant 
angles are not special in this respect (see also appendix A). 

The dependence on the planet mass scale is illustrated in Fig. [4] which shows the evo- 
lution for Q' = 10, twice the minimum masses and without adjustment of the innermost 
semi-major axis. Figure [5] gives the evolution for the same system but with the initial in- 
nermost semimajor axis adjusted to make the innermost period ratio 2.23. Both these runs 
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Figure 1. Time dependent evolution for Q'=10 and minimum masses. Shown is the eccentricity of the innermost planet (upper 
left panel), middle planet (upper right panel) and outermost planet (lower right panel) as a function of time. Running time 
averages of the cosines of the resonant angles are shown in the lower right panel. These correspond to <E>i = 2X2 — Ai — roi(full 
line), <I>2 = 2A2 — Ai — z&2( dotted line), $3 = 2 A3 — A2 — V02 (dashed line), $4 = 2 A3 — A2 — 073 ( dot dashed line). For all plots 
time is measured in yr. 



show similar evolution to the minimum mass case illustrated in Fig. [T] but scaled to larger 
eccentricities. Thus the results are not sensitive to the choice of minimum masses or small 
changes in the semi-major axes. 



3.1 The effect of initial eccentricity 

The calculations described above started with the planets on circular orbits. The motivation 
for this is that it is expected that tidal circularization will lead to a state closely approaching 
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Figure 2. As in Fig[T]but for Q' = 1000, minimum masses and initial innermost period ratio 2.23. 



this provided that initial eccentricities are not so large as to disrupt the system. In order to 
explore this aspect in more detail, we have carried out a number of simulations which started 
with non circular orbits. We have explored initial eccentricities <~ 0.1 and found that the 
system enters into an eccentricity decay mode that causes the system to evolve towards 
the state we found starting from circular orbits in which three of the resonant angles move 
towards libration with non zero values for the running time averages of their cosines. 

In order to illustrate this we describe here results obtained when just one of the planets 
was started with a non zero eccentricity for various values of Q' < 1000. To initiate the system 
it was started as for the simulation illustrated in Figure [T] but with one of the planets set at 
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Figure 3. Evolution of the running time averages of eicos$i (solid line) e2Cos<J>3 (dashed line) and e3Cos<3?4 (dot-dashed 
line) for Q'=1000 and minimum masses. Time is measured in yr. 



pericentre in an orbit with the same angular momentum but with an eccentricity ei mp ~ 0.1. 
This is equivalent to transfer to a less tightly bound orbit with the original specific angular 
momentum, but with the specific energy multiplied by 1 — e? . In all cases the eccentricities 
of the planets become coupled and start to decay. After a while the system attains a slow 
decay mode for which, in a time average sense, the eccentricities are in an approximately 
constant ratio such that e2/e\ ~ 1.2 and es/ei ~ 1.6. 

Fig. [6] illustrates the evolution of the eccentricities and running time averages of the 
cosines of the resonant angles for minimum mass simulations with Q' = 10. Cases for 
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Figure 4. As in Fig[T]but with Q'=10 and with masses a factor of two larger than the minimum. 



which the innermost, middle and outermost planet respectively had an initial eccentric- 
ity are shown. For comparison a case with twice the minimum masses and for which the 
outermost planet was started with a non zero eccentricity is also shown. The evolution of 
all of these systems was as described above. 

However, the introduction of the initial eccentricity while conserving the angular mo- 
mentum of the system causes the semi-major axes and hence the period ratios in the system 
to change. This evolution involves energy dissipation at fixed angular momentum and thus 
the system must spread radially as for an accretion disc (eg. Lynden-Bell & Pringle 1974). 
This causes the period ratios to increase (see Fig. [7]). When this happens the system attains 
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Figure 5. As in Fig[4]but with initial inner period ratio 2.23. 



period ratios that differ from the original ones. However, this is easily compensated for by 
starting with appropriately shifted initial conditions. The important point is that we get 
the same type of evolution independently of such shifts. Fig. [7] shows that the largest shift 
occurs when the outermost planet starts with an initial eccentricity. The period ratio of the 
innermost pair shifts from 2.23 to 2.37. A consequence of this is that the forced eccentricity 
of the innermost planet is ~ 25% smaller in this case than in the case where the inner 
planet started with a non zero eccentricity (see Fig. [6]). We remark that the same evolution 
is obtained when the masses are increased by a factor of two except that the ultimate eccen- 
tricities are a factor between 1.5 and 2 larger. The above trends are in qualitative agreement 
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Figure 6. The evolution of the eccentricities and running time averages of the resonant angles are plotted as functions of 
time measured in years for simulations with Q f = 10. The first second and third rows of panels correspond to the cases where 
the inner, middle and outer planets were started with an initial eccentricity (see text) respectively and minimum masses were 
adopted. The fourth row is the same as the third except that the planet masses were twice the minimum values. From left 
to right for each row, the first panel shows the eccentricity of the innermost planet, the second the eccentricity of the middle 
planet and the third panel shows the eccentricity of the outermost planet. The fourth panel shows running time averages of 
cos($i) (solid line), cos(*3?2) (dotted line) cos(<3>3) (dashed line) and cos(<E>4) (dot-dashed line). 



with a simplified discussion of the resonant interaction between the planets that we give in 
appendix A even though only two of the resonant angles $1 and $3 are found to be librating 
at this point with cos($4) having a non zero time average of relatively small magnitude. 

Finally we have checked that changing Q f gives the same form of solution but with an 
evolutionary time that scales with Q' . Results for Q f = 100 and Q' = 1000 that confirm 
this view when compared with results already described for Q f = 10 are illustrated in Fig. 
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Figure 7. Evolution of the period ratios. The upper curve in each panel corresponds to the period ratio of the inner pair of 
planets and the lower curve to the period ratio of the outer planets. The upper left panel corresponds to minimum masses, 
Q' = 10 and the innermost planet was started with a non zero eccentricity (see text). In the upper right panel the masses were 
twice the minimum values, Q' = 10 and the outermost planet was started with a non zero eccentricity. The lower left panel 
was as for the upper left panel but Q' = 100. The lower right panel was as for the upper right panel but Q' = 100. Time is 
measured in yr. Note that the shifts occurring in the period ratios are independent of Q' . 



[8j However, practical considerations mean that runs with large Q' cannot be followed to 
libration of the resonant angles. Nonetheless when Q' = 100, cos($3) and cos($i) begin 
to attain non zero running means at the end of simulations with minimum masses after 
2.5 x 10 7 yr. This indicates that resonance could be effectively entered for Q' < 10 4 within 
the lifetime of such systems. These results taken together suggest that the planets may be 
undergoing circularization in a three planet resonance. 

As the circularization process involves energy dissipation at fixed angular momentum, 
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Figure 8. The evolution of the eccentricities and running time averages of the cosines of the resonant angles are plotted as 
functions of time measured in years. For the first and second rows of panels, Q' = 100 and Q' = 1000 respectively, minimum 
masses were adopted and the innermost planet was started with a non zero eccentricity (see text). For the fourth Q' = 100 
and the masses were twice the minimum values. In this case the outermost planet was started with a non zero eccentricity. 
From left to right for each row, the first panel shows the eccentricity of the innermost planet, the second the eccentricity of 
the middle planet and the third panel shows the eccentricity of the outermost planet. The fourth panel shows running time 
averages of cos(<J>i) and cos(<I>3) (ultimately uppermost). 



the system must spread increasing the period ratios of neighbouring planets. Although this 
means an increasing departure from strict commensurability, resonant angles remain librat- 
ing or circulating in such a way as to enable the required redistribution of angular momentum 
to take place. A simple model illustrating how this works is given in appendix A. The amount 
of spreading and increase in the period ratios that can take place in a given time depends on 
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Q'. We investigate this in the context of a migration scenario for producing the configuration 
below. 

Simulations incorporating migration indicate that the planets were likely to have been 
in a strict Laplace resonance with all four angles librating. However, if this persists the 
period ratios would increase maintaining a larger departure from commensurability for the 
outer pair than the inner pair in disagreement with observations. Thus some perturbation 
event must have disrupted this resonance enabling the planets to continue circularizing with 
the correct period ratio relationships. As we have seen above, giving one of the planets a 
relatively small eccentricity can produce significant changes to the period ratios. 

To investigate this scenario, we now describe simulations of systems of three planets 
migrating inwards in a disc that set up a Laplace resonance. We go on to study the evolution 
following from its disruption as the result of a perturbation applied once the planets have 
decoupled from the disc. 

4 MIGRATION AND INTERACTION WITH THE STAR: SIMULATION 
RESULTS 

In all the runs we have performed, the planets converge on each other forming a three planet 
4:2:1 resonance, and migrate in locked in that configuration. Whether the resonance can be 
maintained all the way down to the disc inner cavity depends on how fast the migration is. 
There is a tendency for the 4:2:1 resonance to be stable for slow migration rates but to be 
disrupted, with the appearance of higher order resonances once the migration rate becomes 
sufficiently fast. However, the evolution is difficult to predict precisely, as the transition 
from the 4:2:1 resonance to higher order resonances happens through an instability which 
appears to be extended in time and chaotic. Small changes to the model parameters in 
this regime may affect whether this instability occurs or not within a given time interval. 
Similar phenomena have been seen in two planet resonant migration where a sequence of 
resonances can be formed, maintained and subsequently disrupted in turn (see Papaloizou 
& Szuszkiewicz 2009). 

Figure [9] shows the evolution of the semi-major axes of the planets and of the period ratios 
for 3 planets with minimum masses initially located at 0.2, 0.34 and 0.6 au, respectively. The 
disc aspect ratio is H/r = 0.05. Equations ^ and ^ for the migration and circularization 
times are used with Md = 2.5 x 10 -5 M . This was equal to the disc mass within 5 au in the 
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original formalism of Papaloizou & Larwood (2000). However, as there is much uncertainty 
in the migration rates to be used for low mass planets we investigate different possibilities 



using this parameter to scale the migration and circularization rates (see section 2.2). Used 
in this way it no longer corresponds to a disc mass as in the original formulation. 

The initial migration times for this case are then 2.8 x 10 5 and 4 x 10 5 years for the 
innermost and outermost planets respectively. The initial period ratios for the inner and 
outer pairs are 2.21 and 2.34, respectively. This calculation led to the setting up of a Laplace 
resonance that is maintained during the subsequent evolution of the system, although for a 
finite period of a time around 2 x 10 6 years a temporary weak and inconsequential instability 
was present (see below). 

In these calculations , the tidal Q' parameter was set to 10~ 2 . This is at least 3 orders or 
magnitude smaller than the values supposed to apply in Earth-like planets, but this value 
was chosen to give evolutionary timescales that could be reasonably handled numerically 
(i.e. runs that could be done within weeks on a single processor). Tests carried out below 
with more realistic values of Q' indicate that the evolution during the disc migration phase 
is not affected by this specification. Furthermore, when the planets are inside the cavity, the 



scaling given by equation (A37) holds, implying that increasing Q' corresponds to scaling to 
a longer evolutionary timescale. With Q' = 10 -2 , the eccentricity damping timescale given 
by equation (|5| is t s ei = 400 years for a 1 planet at 0.05 au, inducing a timescale on 
which the semi-major axis evolves of about 2 x 10 8 years for an eccentricity of 10~ 3 (see 
section 



2.1). 



Figure [9] shows the evolution of the semi-major axes, the eccentricities, the period ratios 
and the resonant angles $i = 2A 2 — Ai — £j\ and $ 3 = 2A 3 — A 2 — u) 2 . As pointed out above, the 
Laplace resonance for which all four resonant angles librate is established while the planets 
migrate in the disc, after about 5 x 10 5 years. The resonant angles $i and $3 stay close to 
zero. Then, while some of the planets are still in the disc, a weak instability appears for a 
finite period of time (after about 2 x 10 6 years) and the resonance is weakened with the 
angles driven to large amplitude libration with an occasional circulation. But the instability 
eventually dies out and the resonance is recovered before the planets decouple from the disc. 
Note that this instability is not always seen in the runs we have performed. As indicated 
above it seems to be associated with an approach to a transition to a resonance of higher 
degree. This is indicated by the fact that when the migration rate was speeded up by a 
factor of 4, the outer pair of planets underwent a transition to a 3 : 2 resonance while the 
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Figure 9. Semi-major axis a (in au, plots a and 6), eccentricity e {plot c), period ratios (plot d), resonant angles <I>i = 
2\'2 — Ai — oil (in rad., plot e) and $3 = 2A3 — A2 — U12 (in rad., plot f) versus time (in million years) for 3 planets initially 
located at 0.2, 0.34 and 0.6 au. On plots a and b, each solid line represents a planet, whereas the dotted line indicates the 
location of the disc inner edge. Note the different time interval for plot b. On plot c, the solid and dotted lines represent the 
period ratio of the inner and outer pair of planets, respectively. Here = 2.5 X 10 -5 Mq and Q' = (note that the 

evolution during the disc migration phase is not sensitive to the value of this parameter). Even though strict commcnsurability 
is lost as a result of tidal interaction with the star, the planets stay locked in resonances, as indicated by the behaviour of the 
resonant angles. 



inner pair maintained a 2 : 1 commensurability at the point at which the semi-major axes 
had attained the same configuration. 

During the migration phase, the eccentricities are pumped up to several hundredths. 
Because of the low value of Q', tidal interaction with the star begins to affect the evolu- 
tion of the system slightly before the outermost planet penetrates inside the cavity. Once 
the outermost planet reaches the disc inner edge, the disc is removed altogether from the 
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calculation, so that if the planet is pushed out again it will not interact with the disc. This 
procedure is introduced to take into account the disappearance of the disc on a timescale 
of a few million years. In this particular run, the outermost planet enters the cavity after 
3.5 x 10 6 years. At that point, changes in semi-major axes are due only to tidal interaction 
with the star. As the eccentricities of the two inner planets are damped from a few 1CT 3 to 
1CT 3 , their semi-major axes decrease by about 1 percent (see plot b of fig. [9]) in such a way 
that their period ratio increases by somewhat more than a percent, going from 2 to 2.037 
(see plot d of fig. [9]) . 

In the meantime, the eccentricity of the outermost planet is hardly affected, as it was 
already very small (a few xlO -4 ) when the planet entered the cavity. Because of the con- 
servation of the total angular momentum though, its semi-major axis increases by about 
1 percent, which results in the period ratio of the outer pair of planets increasing by a few 
percent, from 2 to 2.077. Therefore, strict commensurability is lost. However, we note that 
the resonant angles still librate, so that the Laplace resonance is not disrupted. The fact 
that the deviation of the outer period ratio from 2 is twice the deviation of the inner period 



ratio from 2 is a consequence of the Laplace relation (equation A28 of appendix A) and 



this feature is preserved in the subsequent evolution of this and similar cases. As we shall 



see, these deviations increase approximately oc t 1 / 3 , as predicted from equation (A37). 

As noted above, the Laplace relation is not satisfied in HD 40307. As we shall see below, it 
is possible to obtain the correct relationship between the period ratios provided the Laplace 
resonance is disrupted such that only $i,$3 and $4 librate and contribute to the secular 



evolution (see also section 3.1) . 

In order to investigate the dependence on the specification of the circularization rates 
induced by the disc, we have performed simulations for which these depend only on the 
planet mass and are therefore constant. We here consider examples for which: 



t m!i = 7 x 10 6 yr, 8 
rrii 

and: 

<i = 8 x 10 3 ^ yr . (9) 
Results for cases both with the minimum masses for HD 40307 and twice the mimimum 



masses are plotted in Figure 10 During the disc induced migration, a three planet resonance 



is formed in both cases. As expected, the evolution for the higher mass case is twice as fast 
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but otherwise similar. The resonance shows only small signs of instability while the planets 
migrate inwards. This is probably helped by the fact that the migration rate is here constant 
rather than increasing inwards as in the previous simulation. As the migration continues, 
the resonant angles $j (i = 1, ...,4) all enter libration, but after the planets enter the cavity 
there is an adjustment so that only $1 and $3 librate. The system then evolves back towards 
a Laplace resonance where all angles librate or contribute to long term time averages. This 
latter aspect is similar to what was seen in the simulations with larger Q' illustrated in 
Figs.[TH5} 

This final evolution is one where the eccentricities all reduce and the planets physically 
separate with increasing period ratios for the inner and outer pairs of planets. This form of 
evolution is to be expected in a Keplerian system in which energy is dissipated while the 
total angular momentum is conserved, as is the case for an accretion disc (e.g., Lynden-Bell 
and Pringle 1974). It is this form of evolution we investigate below and how far the tidal 
effects of circularization can bring the system towards the HD 40307 system. 

4.1 Disruption of the resonance 

As the continuation of the simulations described above maintains a three planet Laplace 
resonance, the period ratios will not have the same form as in HD 40307. As indicated 
above the system has a larger period ratio for the inner pair as compared to the outer pair 
whereas the Laplace relation leads to the converse. However, we have found that there are 
resonant interactions for which only $!,$ 3 , and $4 librate. Hereafter we refer to this as 
the standard state of libration. This can lead to the required form for the period ratios but 
some disruption of the form of the resonance initially set up within the cavity needs to take 
place. Mechanisms could involve close encounters or collisions occurring shortly after the 
disc migration phase. 

We comment that we have found it impossible to disrupt the 2:1 resonance of the inner 
pair of planets with a slow process. For instance, we have added a fourth planet on an inclined 
orbit, but found that even a massive perturber (3 Jupiter masses) at 5 au could not take 
the inner pair away from the 2:1 resonance. To do this, the parameters of one or more of the 
orbits have to be varied suddenly preventing the resonance from responding adiabatically. 
An encounter with a fourth planet on a parabolic orbit is a possible means to achieve this and 
there are other possibilities involving one or more direct collisions. To affect the resonance, 
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Figure 10. The evolution of model systems with the minimum masses of HD 40307 (left panels) and twice the minimum 
masses (right panels) showing the formation and maintenance of three planet resonance as a result of disc migration. For these 
models, the migration and circularization rates are constant and Q' = 10 (although evolution during the disc migration phase 
is not sensitive to this parameter) 



The upper panels give the evolution of the semi-major axes and the lower panels the period ratios. 



the semi-major axis of one of the planets should change by more than an effective libration 
width associated with one of the resonant angles. We have performed simulations for which 
the velocity components of one of the planets were all suddenly increased or decreased by a 
constant factor. This has the effect of preserving circular orbits which is reasonable under 
conditions of effective circularization. We found that fractional changes exceeding ~ 2 x 1CT 3 
were sufficient and, provided the changes were not too excessive, simulations approach the 
same form of evolution. 
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Below we describe a simulation where an impulsive change that decreased the velocity 
components of the innermost planet by 0.996 was applied just after the outermost planet 



entered the cavity. The end of the simulation illustrated in Fig. 10 (or its larger planet mass 
counterpart where appropriate) was used as a starting point for the simulations presented 
below, as these are reasonable outcomes for a system having undergone inward migration 
induced by a disc and then entered an inner cavity. In some cases, the tidal parameter Q' 
was also adjusted. 

We remark that we have performed simulations with other types of applied impulse. 
For example we have augmented the above form by giving the innermost planet an initial 



eccentricity in a similar manner to that undertaken in section 3.1 When this is small and 
typically of order 10~ 3 this makes no difference on account of rapid effective circularization. 
We have also considered impulsive changes to the outermost planet giving it an eccentricity 
~ 0.1 as in section 3.1| In this case period ratios can be shifted to be close to the observed 
ones just after the outermost planet has entered the cavity (see below). 



4.2 Evolution of systems on three planet resonances undergoing 
circularization and radial spreading 



In Fig. [TTJ we illustrate the evolution of a model system of three planets with masses twice 
the minimum masses of HD 40307 with Q' = 10. Two runs are displayed. The first one 
continues the evolution of the corresponding case plotted in Fig. [10] In the second run, an 
impulse was first applied to the innermost planet that reduced all velocities by a factor of 
0.996. 



As shown in Fig. [TTJ during the subsequent evolution the system spreads radially with 
increasing period ratios. The way this happens is illustrated by the simple model described 
in appendix A. In this discussion the period ratios are replaced by the quantities = 
n i/ n 2 — 2 and = n 2 /n^ — 2 that are easily obtained from Fig. 11 We have been 



able to find good fits to the time evolution of the form found in appendix A as given by 



equation (A37) 



a 



1,2, 



(10) 



where a and /3 are constants and t$ is the time in units of 10 6 yr measured from the start 
of the relevant simulation shown in Fig. 10 

In the case with no initial impulse, a = 2.77 x 10 5 and (3 = 1.5 for the outer pair, 
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Figure 11. Evolution of a model system of three planets with masses twice the minimum masses of HD 40307. The left panels 
correspond to a continuation of the evolution from the corresponding case illustrated in Fig. |10| The right panels correspond 
to the evolution after an impulse was applied to the inner planet (see text). For both models Q' = 10. The period ratios ni/712 
and n2/ri3 are plotted as a function of time in the upper panels. The upper curve corresponds to the outer pair (left panel) 
and the inner pair (right panel). In this and similar figures below, the triangles correspond to analytic fits of the form given 
by equation 1 10| (and see text). The lower panels show the evolution of the eccentricities of the three planets. In the lower left 
panel the curves from top down correspond to the inner, middle and outer planets respectively. In the lower right panel the 
curves from top down correspond to the middle, lower and outer planet. 



and a = 2.37 x 10 6 and =1.5 for the inner pair. For the case with an initial impulse, 
a = 7.73 x 10 6 and /3 = —26.5 for the outer pair, and a = 9.66 x 10 5 and (3 = —26.5 for 
the inner pair. For the latter case, a standard libration state occurred with the deviation of 
the inner period ratio from two being about a factor of two larger than the corresponding 
deviation for the outer period ratio, similar to the situation in HD 40307. This is in contrast 
to the case with no impulse that persisted in a Laplace resonance with all angles librating. 
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In that case, the Laplace relation persists which leads to the relative deviation of the period 
ratios from two behaving differently from HD 40307 (see above). 

Despite this, we can estimate the time it would take to obtain the outer period ratio of 
2.13 appropriate to HD 40307. This is 6.1 x 10 8 yr. The scaling of the circularization time, 

2/3 

t s ei oc m i leading to the scaling of the evolution time, t, to reach a given state, with mass 
of the form t oc m i 8//3 , implies that for a minimum mass system this would be ~ 4 x 10 9 yr. 
Thus, it might just be possible to achieve the outer period ratio in this case but then the 
inner period ratio would have to be determined subsequently by alternative means. 

For the case with the applied impulse, we may estimate the time to bring the inner 
period ratio to the appropriate value 2.23. This is somewhat uncertain because of the small 
variation seen in the simulation, but it is about ~ 1.1 x 10 10 yr. Note that the relative 
period ratios are not precisely identical to those of HD 40307. In particular the ratio of the 
deviation of the inner period ratio from two and, the deviation of the outer period ratio from 
two is about 30% too large. However, this may be adjusted by changing simulation input 
parameters. For example we recall that the simulations reported here were performed with 
the same value of Q' for all of the planets. To indicate one possibility we note that tests we 
have performed indicate that the above ratio may be reduced by taking a smaller value of 
Q' for the outermost planet than the others. 

We remark that, while it might be barely possible to obtain an inner period ratio of 2.23 
for a larger mass system, it is almost certainly not possible for the minimum mass system 
as the time required would be ~ 6.98 x 10 10 yr. 

Note that, as the evolution progresses, the eccentricities of all the planets tend to de- 
crease, as does their variation. This is a general consequence of the spreading that is driven 
by the energy dissipation in these simulations. 

To investigate further the form of the evolution and its scaling with run parameters, we 



plot in Fig. [12] the results of two simulations of models with minimum masses and for which 
an impulse of the above form was applied to the innermost planet. We compare cases with 
Q' = 0.01 and Q' = 0.1. As stated above, we had to consider rather small values in order 
to perform simulations that produced a reasonable variation in a reasonable time. These 
resulted in a standard state of libration and relative period ratios very similar to the higher 



mass simulation with an initial impulse illustrated in Fig. 11 



In the case with Q' = 0.01, the fit parameters are a = 3.58 x 10 3 and /3 = 1.7 for the inner 
pair, and a = 4.20 x 10 4 and =1.7 for the outer pair. For the case with Q' = 0.1, they 
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Figure 12. As in Fig. |11| but for models with minimum masses and for which an impulse was applied to the innermost planet 
(see text). For the left panels, Q' = 0.01 and for the right panels, Q' = 0.1. In the upper panels, the upper curves correspond 
to the inner period ratio. In the lower panels, the curves from top down correspond to the middle, inner and outer planets 
respectively. 



are a = 4.13 x 10 4 and (5 = 1.5. for the inner pair, and a = 4.18 x 10 5 and (3 =1.5 for 
the outer pair. The values of a, being about an order of magnitude larger in the second 
case, are consistent with the evolution time being an order of magnitude larger, thus scaling 
with Q' as expected. The time to attain an inner period ratio of 2.23 in the Q' = 0.01 
case is ~ 4.4 x 10 7 yr, leading to an estimate of ~ 4.4 x 10 10 yr when Q' = 10. This is 
less than the estimate based on the analytic scaling with planet mass implied by equation 



(A37) applied to the higher mass simulations illustrated in Fig. 11 by a factor ~ 1.5, which 
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we believe is a reasonable level of consistency in view of the extent of the extrapolation 
involved. The eccentricities in these simulations decrease with time in a manner analogous 



to the simulations shown in Fig. 11 



The long times required for the above cases to reach a configuration similar to that now 
observed in HD 40307 arise because only a very small perturbation was applied just after 
the systems decoupled from the disc, requiring a large amount of subsequent period ratio 
evolution. Such a scenario is only feasible for masses exceeding twice the minimum and small 
Q' ~ 10. 

However, if a larger perturbation is allowed, it can provide a larger shift in the period 
ratios requiring less evolution later on (see Fig. [7]). This is equivalent to starting on the 



evolutionary curves shown in Figs. [Tlp^ at a later time. 

71 was 



For example we found that if a perturbation of the type adopted in section |3.1 
applied to the outer planet a shift of the right type could be obtained. We found that when 
the imposed eccentricity of the outer planet was taken to be 0.15, the period ratios were 
shifted to very similar values to those in HD 40307. Thus if large enough perturbations are 
allowed, only a small amount of subsequent evolution may be needed. 



5 DISCUSSION AND CONCLUSION 

In this paper, we have shown that the planets around the star HD 40307 could be undergoing 
a resonant interaction, despite departure of the period ratios from very precise commensu- 
rability. This is indicated by the fact that three of the four resonant angles librate or are 
associated with long term changes to the orbital elements. Note that such a resonant state 
can occur without exact commensurability only because the eccentricities are very small, 
which can be brought about as a result of tidal circularization from a state with initially 



higher eccentricities (see section 3.1). When this occurs in a system that starts from a near 
2:1 commensurability, the difference of the period ratios from two then increases oc 

We propose that the planets in this system were in a strict Laplace resonance while they 
migrated through the disc, with all 4 resonant angles librating. Exact commensurability was 
then departed from as indicated above as a result of the tidal interaction with the star, 
which preserved libration of at least some of the resonant angles, after the planets entered 
the disc inner cavity. Because of the Laplace relation , the period ratios evolve in such a way 
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that the ratio for the inner pair of planets is smaller than that of the outer pair of planets. 
The opposite is observed in HD 40307. 

To get the appropriate form for the period ratios, some disruption of the resonance set 
up in the cavity needs to take place (this is in contrast to claims by Zhou 2009). A close 
encounter for instance might be a possible cause. An impulsive interaction would be sudden 
enough that the resonance would not respond adiabatically. If the velocity of the innermost 
planet is reduced by a factor 0.996 for instance, deviation of the inner period ratio from two 
gets twice as large as the corresponding deviation for the outer period ratio, while 3 of the 
4 resonant angles still librate, similar to the situation in HD 40307, within the lifetime of 
the system. 

However, we found that we could get a period ratio of 2.23 for the inner pair only if the 
masses are significantly larger than the minimum masses and tidal dissipation is possibly 
unrealistically effective. From the results obtained in section |4j when the Laplace resonance 
is disrupted at an early stage, we estimate that the deviation of the inner period ratio from 
two is approximately given as a function of time by W 1>2 ~ 0.23(0. 157mi /3 t 1 o/Q , 10 ) 1/3 , where 
Qw = Q'/10> *io — t/(10 10 yr) and mi = mi/4.2 M e (m 2 /mi and m 3 /mi are assumed fixed 
for the purposes of this discussion). Accordingly, to attain W12 = 0.23, we would require 
masses exceeding the minimum estimate by about a factor of two for Q' 10 ~ 1. 

If the Laplace resonance is preserved, from the results obtained in section |4j we may also 
estimate that the deviation of the outer period ratio from two is approximately given as a 
function of time by ~ (2Mm^ 3 ti /Q' l0 ) 1 ^ 3 . Thus it may be possible to get a period ratio 
of 2.13 for the outer pair of planets with near to minimum masses in a reasonable time if 
Q' l0 ~ 1 for all the planets, but then the inner period ratio would be only ~ 2.065. It would 
thus seem likely that additional stronger dynamical interactions are required to bring the 
system to its final state. 

For example we found that if the outer planet was given an eccentricity e% ~ 0.15 shortly 
after decoupling from the disc, period ratios similar to those observed in HD 40307 could be 
produced on a time scale of a few x 10 6 yr. (see sections 
in this paper indicates that tidal circularization would continue through resonant interaction 
with the period ratios separating after these interactions are over. 

Note that a similar type of scenario may apply to the system of planets around GJ 581, 
which is currently closer to exact commensurablilty than HD 40307. Indeed, different mean 
motion resonances can be established while the planets migrate within the disc, depending 



3.1 and p| . Then, the work presented 
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on the initial relative location of the planets and on the migration timescale. However, as 
it appears that this system may have been in higher order resonances than HD 40307, it 
is possible that these would have been broken completely as a result of circularization. We 
intend to present a study of GJ 581 in a forthcoming paper. 

In a study of the tidal evolution of the system around HD 40307, Barnes et al. (2009) 
have concluded that the planets cannot be Earth-like. According to them, if it were the 
case, extrapolating back in time the tidal evolution of the eccentricities from current values 
(estimated from direct simulations of the observed system), the system would have been 
unstable. The work in this paper does not support such a conclusion. We note that Barnes 
et al. (2009) neglected the dynamical interactions between the planets, which clearly will 
have played a role if the current eccentricities are small because resonant interactions coupled 
with tidal effects have been associated with long term changes to the orbital elements. 

The results reported in this paper indicate that resonances in multiple low mass-planet 
systems may play more of a role than currently claimed in the literature. Significant de- 
parture from exact mean motion resonance does not necessarily preclude a resonant state 
when, for example, tidal effects are important causing eccentricities to be small, as libration 
of some of the resonant angles in that case may still exist, and be associated with long term 
evolution of the orbital elements. Although the confirmation of resonance effects is prob- 
lematic for small eccentricities, we expect that further detections and analysis of systems 
similar to HD 40307 will test the relevance of the scenarios explored in this paper. 

Finally, we comment that, because of the rapid increase of the effectiveness of tidal 
interactions with decreasing orbital period, the dynamical interactions discussed here would 
be of much greater importance for systems in which the period of the innermost planet is 
significantly shorter than in the HD 40307 system. For example, we note that the planet 
CoRoT-7 b, which has a projected mass of 11.12 M®, has a period of only 0.85 day (Leger 
et al. 2009). If this short period orbit is a result of migration, it is likely this planet was 
pushed inwards by a companion which would have been on a commensurable orbit (Terquem 
& Papaloizou 2007). We note in support of this idea that pre-main sequence stars typically 
have rotation periods of several days (Bouvier et al. 1993), which would imply that the disc 
inner cavity should extend well beyond the orbit of CoRoT-7 b if that were magnetically 
maintained. In that case, this planet could not have migrated so far inwards without being 
shepherded by a companion. 
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APPENDIX A: SEMI-ANALYTIC MODEL FOR A SYSTEM IN A THREE 
PLANET RESONANCE UNDERGOING CIRCULARIZATION 

In this appendix we develop a semi-analytic model that shows how a two or three planet sys- 
tem undergoes radial spreading as a result of evolution driven by tidal circularization. The 
coupling between the planets occurs because resonant angles may librate or be associated 
with long term changes even though there may be significant deviations from exact commen- 
surability. We begin with the system without tides which is Hamiltonian. Subsequently, we 
add the effects of tides. 



Al Coordinate system 

We adopt Jacobi coordinates (Sinclair 1975, Papaloizou & Szuszkiewicz 2005) for which 
the radius vector of planet i, r*, is measured relative to the centre of mass of the system 
comprised of M and all other planets interior to i, for % = 1,2..., N with N = 3 here and 
% — 1 corresponding to the innermost planet. The Hamiltonian, correct to second order in 
the planetary masses, can be written in the form: 

3 



H — > -m, • r 
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Here Mj = M + rrii and = r j — r j . 

The equations of motion for motion in a fixed plane, about a dominant central mass, 
may be written in the form (see, e.g., Papaloizou 2003, Papaloizou & Szuszkiewicz 2005): 



dH 

E, = (A2) 



i dH dH 

= aI" + "-aB- (A4) 

*• = 3V (A5) 
Here the orbital angular momentum of planet i which has reduced mass mi = m^Mj (M + 
rriio), with rriio being the actual mass, is Lj and the orbital energy is E'j. For motion around 



32 Papaloizou and Terquem 
a central point mass M we have: 



Li = m^GMiaiil-ef), (A6) 
E, = -<^, (AT) 

where = M + m^o, a« denotes the semi- major axis and e« the eccentricity of planet i. 
The mean longitude of planet i is Aj = rij(£ — £ 0i ) + where rij = a/ GMija\ is its mean 
motion, with £ 0i denoting its time of periastron passage and W; L the longitude of periastron. 
We remark that results valid for a two planet system, i = 1,2, may be obtained by setting 

f«3 = 0. 

The Hamiltonian may quite generally be expanded in a Fourier series involving linear com- 
binations of the five angular differences w\ — w 2 , ru 2 — 073 and Aj — zui, i — 1,2,3. 

Here we are interested in the effects of two first order 2 : 1 commensurabilities associated 
with the outer and inner pairs of planets. In this situation, we expect that any of the four 
angles $1 = 2A2 — Ai — w\, $2 = 2A2 — Ai — w 2 , $3 = 2A3 — A2 — 'cu 2 and $4 = 2A3 — A2 — s7 3 
may be slowly varying. Following standard practice (see, e.g., Papaloizou 2003, Papaloizou 
& Szuszkiewicz 2005), only terms in the Fourier expansion involving linear combinations of 
$1, $ 2) $3 and $ 4 as argument are retained. Working in the limit of small eccentricities 
applicable here, only terms that are first order in the eccentricities need to be retained. 
Expanding to first order in the eccentricities, the Hamiltonian may be written in the form: 



where: 



H = E 1 + E 2 + E 3 + H 12 + # 23 , (A8) 



Gun ■ TYi ■ 

= — [ejCij cos($ i+i _i) - e<Aj cos($ 2 i_i)] , (A9) 
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+ Zb\ /2 (x iJ )-4x iJ , (AlO) 



+ 46? /2 (x y ) I . (All) 



Here b^^ 2 (x) denotes the usual Laplace coefficient (e.g., Brouwer & Clemence 1961) with the 
argument Xij = ai/aj. 
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A2 Behaviour near a resonance with disc tides incorporated 



Using equations (A2)-(A5) together with equation (A8) we may first obtain the equations 
of motion without the effect of circularization tides and then add this in. We obtain correct 
to lowest order in the eccentricities: 
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(A12) 
(A13) 
(A14) 
(A15) 



(A16) 
(A17) 

(A18) 
(A19) 
(A20) 
(A21) 



0,1,0,2,0,3 



The terms on the right hand sides involving the circularization times t s e i describe the effects 
arising from tides raised by the central mass on planet i. The terms oc ef/t s ei in equa- 



tions (A15), (A17) and (A18) account for the orbital energy dissipation occurring as a result 
of circularization at the lowest order in ei for which it appears. 



A3 Energy and angular momentum conservation 

In the absence of tides, the total energy E = H and angular momentum L = L\ + L 2 + L 3 
are conserved. When circularizing tides act, the total angular momentum is conserved but 
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energy is lost according to: 

dE (A22) 



dt f-f t*. 

i=i e ' J 



A4 Resonance tides and migration 

When both pairs of planets (1, 2) and (2, 3) are in resonance, a three planet resonance is said 



to exist. It is convenient to perform a time average of equations (A12)-(A21) taken over a 
time long compared to both the characteristic orbital period and the libration period , but 
short enough compared to the tidal time scale (see eg. Sinclair 1975 for a discussion of this 
aspect). The angles <3?j may be undergoing libration or circulation. 



Equations (A12)-(A17) may be averaged directly, while equations (A18)-(A21) may 
be averaged directly for small amplitude librations or alternatively after multiplying them 
by quantities such as either, e.g., cos$i, cos $2, cos $3 and cos $4 or eicos$i, e 2 cos$ 2 , 
e 2 cos $3 and cos $4, respectively. The time averages of these quantities as well as sin $j 
and ejsin$j, (i = 1,2,3) are assumed to vary on the tidal time scale or longer. When 
these are non zero, long term changes to the orbital elements dj and may occur through 
planet-planet interactions even without tides. As illustrated in Figures [2] and [3j six of these 
quantities are found to exhibit steady long term time averages of the required type that are 
quickly established for a model of HD 40307 with assumed Q' = 1000. 

As this averaging procedure results in equations of essentially the same form as the case 
where the participating resonant angles (defined to be those that are associated with non 
zero time averages) are assumed to have averages close to but not equal to zero or tt, while 
undergoing small amplitude librations, leading to the same scaling with physical parameters, 
but with an adjustment of numerical coefficients. For simplicity we restrict consideration to 



that case. Assuming small amplitude librations of all four angles, equations (A18)-(A21) 
give relations of the form: 

2n 2 -n 1 = — , (A23) 

eiMia 2 

2n 2 -m = /3, (A24) 

2n 3 -n 2 = 0, (A25) 
m 2 n 3 C 23 cos $ 4 

2n 3 -n 2 = — , A26 

e 3 M 3 

where: 

_ m 3 n 2 a 2 D 23 cos $ 3 _ m x n 2 C 12 cos $ 2 

e 2 a 3 M 2 e 2 M 2 1 ' 



Multiple systems of super-Earths 35 
and the in the above can be set to either zero or n. 



We note that, when all four angles participate, equations (A24) and (A25) imply the 
strict Laplace resonance relation 

3n 2 - 2n 3 - n x = 0. (A28) 
However, simulations of HD 40307 indicate that only $i, $ 3 , and $4 participate. Thus equa- 



tion (A24) is absent and the strict Laplace relation does not hold. Indeed, this relation 
implies that C = (rii/n 2 — 2)/(n 2 /n 3 — 2)n 2 /n 3 should be unity. Being about 4, the Laplace 
relation is not satisfied. When $ 2 does not participate and the Laplace relation does not 
hold, the quantity C is not constrained from the outset but is determined as a consequence 
of the initial conditions and following time dependent evolution. 

Below, we consider the case when either all four angles contribute or only $1, $ 3 and $ 4 
contribute. We note that a similar discussion applies when only $i,$ 2 and $ 3 participate 
that leads to the same conclusions (e 3 decays independently due to tides in that case and 
may be taken to be zero). In all cases, one can readily verify that the above system provides 
a complete set of equations for the time averages of the resonant angles and the orbital 
elements when the form of the tidal forces is specified. 

However, we note that we can obtain the characteristic form and time scale of the evo- 



lution by considering the conservation energy given by (A22) together with the constancy 



of the total angular momentum. To do this we begin, by noting that equations (A23) and 



(A26) provide equations for the mean motion ratios of the form: 
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(A29) 



(A30) 



and /2,i and / 3j2 are defined implicitly through equations (A23) and (A26). In addition, 



equation (A25) gives e 2 in terms of e 3 , £2,1 and £ 3)2 in the form: 

e 2 m 3 M 3 n2a2D 2; 3 cos $ 3 — miM 3 n 2 a3Ci t 2 cos $ 2 
e 3 m 2 M2n 3 a 3 C2 3 cos $ 4 



(A31) 



When cos$ 2 has zero time average and so is replaced by zero, equation (A24) does not 
hold and accordingly the Laplace relation also does not hold. We comment that £ 2j i and 
£ 3j2 measure the deviation from strict commensurability of the corresponding mean motions 
and do not have to be small. For a system such as HD 40307 these are large, particularly in 
comparison comparison to the characteristic size of the eccentricities. 
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We remark that according to equations (A23) - (A26), other things remaining fixed, the 
eccentricities are proportional to the planet mass scale and, inversely proportional to the 



deviation from commensurability. The results presented in section |3.1| confirm these trends 
but with some deviation that is to be expected because all the relevant resonant angles, 
although contributing to the evolution, are not in a state of small amplitude libration in 
those simulations. 



A5 Energy and Angular momentum conservation 

Replacing Mj (i = 1, 2, 3) by M to simplify expressions without loss of content (as M ^> raj), 
the total energy and angular momentum may be written, respectively, as: 
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(A33) 



mi /J( 3 mi(/ 3 , 2 / 2 ,i) 1/3 _ 
Here, we have neglected in L terms proportional to the squares of the eccentricities as 
these can be shown to lead to small effects compared to those arising from the eccentricity 



dependence of / 2 ,i and / 3j2 . We now use the above together with equation (A22) to form the 



quantity L 2 d(EL 2 )/dt = dE/dt + (2E/L)dL/dt and find that this is given by: 
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We may write equation flA34ft in terms of the deviations from commensurabiity: W 2 1 
ni/n 2 — 2 = f 21 — 2 and W 2j3 = n 2 /n 3 — 2 = f 32 — 2, obtaining: 
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In addition, the right hand side can be expressed in terms of Wi >2 and W 2>3 using equa- 
tions ( |A23D ( |A26| ) which give: 



E 



i=l 



2e 2 Ei Gmim 2 l n\a 1 D'l 2 cos 2 $i 



1 + 



e 2 \ 2 -E 2 te,3 



Gm 3 m 2 n 3 C 2 3 cos $ 4 



(A36) 



e,3 



e 3 / -c/3t e , 2 

We may now obtain characteristic form of the time evolution of the departure from strict 
commensurability. To do this, we may simply assume that W 2> x, W 2}3 are comparable, signifi- 
cantly less than unity and ~ Wo, say, as has been found numerically in general and must be at 
least approximately the case for a strict Laplace resonance for which W 2j s = 2W 2j i/(l — W 2t i) 
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(see equation A28| ). In addition, we take the planet masses to be comparable and 



m. 



Then the characteristic evolution expected from equations (A35) and (A36) can be obtained 
from an equation of the form dWo/dt ~ (m/M) 2 / (t c Wo ), where we take t c to be the shortest 
circularization time t P ,-. Thus: 
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(A37) 



where to specifies an arbitrary origin of time. We recall that (A37) applies to a two planet 
system by appropriately setting 777,3 = as well as a three planet system. The time required 
to attain a departure Wq is t — to ~ 3Wq (M/m) 2 t c . In general, this is very much longer than 
the characteristic circularization time. 
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